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As a complement to experimental and theoretical approaches, numerical modeling has become an 
important component to study asteroid collisions and impact processes. In the last decade, there have 
been significant advances in both computational resources and numerical methods. We discuss the present 
state-of-the-art numerical methods and material models used in ’’shock physics codes” to simulate impacts 
and collisions and give some examples of those codes. Finally, recent modeling studies are presented, 
focussing on the effects of various material properties and target structures on the outcome of a collision. 


1. INTRODUCTION 

The modeling of impact processes can be based upon 
mathematical synthesis of experimental results, on direct 
theoretical application of the principles of physics, or on 
the use of those principles in numerical codes. 

The direct application of experimental results is not usu¬ 
ally possible, because the experiments cannot be performed 
at the actual conditions of interest. To bridge the gap, scal¬ 
ing theories are developed using physical principles to ex¬ 
trapolate experimental results to the actual conditions of in¬ 
terest. For some time the principal scaling theories have 
been based upon the physical concept of a ’’point source”, 
wherein the genesis of an impact process is considered to 
occur instantaneously at a negligibly small region on the 
surface of the target object. Prior to 1982 the point source 
scaling was assumed to be governed by the kinetic energy 
of the impactor, but Holsapple and Schmidt (1982) showed 
that such an assumption was not warranted, and extended 
the analysis to arbitrary point sources. Their general ap¬ 
proach has been followed in numerous subsequent papers 
by Holsapple, Housen, and Schmidt, one can refer to the re¬ 
view of those scaling approaches in Holsapple (1993). That 
scaling theory continues with applications to date. 

In principle, the physics that governs such processes is 
known, at least in the continuum approach. That physics 
includes the balance laws of mass, momentum, and en¬ 
ergy, augmented by mathematical descriptions of the ma¬ 
terial behavior. Material behavior commonly includes the 
’’equation of state”, which models the hydrostatic compo¬ 
nents of the stress histories and the principal thermodynam¬ 
ics; as well as equations describing the deviatoric (shear) 


components of the response. The latter descriptions include 
stress-strain-temperature-relations and include also models 
for failure, flow, and fracture. That material behavior is the 
source of the primary uncertainties about the correct way 
to model these processes. However, there has been much 
progress in the last couple of decades, so that direct numer¬ 
ical solutions using time and space-stepping increments are 
becoming increasingly sophisticated and important. 

The codes used in those numerical approaches are com¬ 
monly called ’’hydrocodes”, a remnant of the early days 
when they were used in the military industry to make cal¬ 
culations of weapons effects, and the modeling only in¬ 
cluded the hydrodynamic aspects of the processes. Nowa¬ 
days those codes are better called ’’shock-physics codes” 
(Pierazzo et ai, 2008). 

Just within the last decade, another approach has been 
applied to asteroid processes. That approach has been bor¬ 
rowed from the fields of particle mechanics and of ”n-body” 
studies. They model the material of the body as a large 
number of individual, discrete, usually spherical, usually 
mono-sized, indestructible particles. With those simplifica¬ 
tions the balance of momentum alone determines the mo¬ 
tions of the particles. Balance of mass is automatic, and 
there is no accounting for energy balance and heating. The 
interactions of the particles is modeled using combinations 
of concepts of restitution, friction, and most recently cohe¬ 
sive forces, with a number of interaction parameters. These 
approaches have only been made possible because of the 
extensive growth of computing power. However, because 
of the extreme complications of the interactions of real par¬ 
ticles at high energies and stresses, those approaches will 
most likely be restricted to cases where the stresses remain 
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low, and the discrete particle nature of the process can iden¬ 
tify processes not contained in the continuum approaches. 

In this chapter, we give an overview of the important 
asteroid properties that determine the outcome of a colli¬ 
sion and discuss the physical processes involved. We then 
present recent experimental results and the theoretical ap¬ 
proaches to describe the outcome of an impact. The main 
part of this chapter is devoted to a detailed discussion of the 
current state-of-the-art numerical models - shock physics 
codes - which are used in the field to simulate impacts 
and collisions. Some examples of those codes are pre¬ 
sented and the various approaches to model the important 
properties and processes are detailed here, including hybrid 
hydrocode-particle code computations. Finally, some ex¬ 
amples of recent modeling are presented. In these studies, 
the effects of various material properties and target struc¬ 
tures on the outcome of a collision are discussed. 

2. IMPORTANT PROPERTIES AND PROCESSES 

Asteroids have complex shapes, internal structures and 
material properties. The impact response and mechanical 
behavior of such objects is naturally difficult to model. In 
this section, we present some of the important internal aster¬ 
oid properties that determine the outcome of a collision and 
discuss late-stage processes. Asteroid interiors, morpholo¬ 
gies and surface geophysics are further discussed in Scheers 
et al., Murdoch et al. and Marchi et al. (this volume). 

2.1 Porosity 

The fact that the bulk density of many asteroids is well 
below the grain density of their likely meteorite analogues 
indicates that many have significant porosity ( Britt et al., 
2002). In particular, several lines of evidence point to 
the presence of a high degree of porosity for asteroids 
belonging to the C taxonomic class, as indicated by the 
very low bulk density (~ 1.3 g/cm 3 ) estimated for some 
of them, such as the Asteroid 253 Mathilde encountered by 
the NEAR Shoemaker spacecraft ( Yeomans et al., 1997), 
and as inferred from meteorite analysis ( Britt et al., 2006). 
That porosity might be a result of a rubble-pile structure, 
as suggested for instance by spacecraft observations of as¬ 
teroid Itokawa ( Fujiwara et al. 2006). The various forms 
of structures and porosities of asteroids were discussed in 
Asteroids III ( Britt et al., 2002; Richardson et al., 2002; 
Asphaug et al., 2002). 

It is useful to distinguish between ’’micro” and ’’macro¬ 
porosity”. The distinction is primarily a matter of scale. 
The terminology arose in the study of meteorites, wherein 
porosity not apparent to the naked eye was called micro¬ 
porosity, and visually obvious voids between a grain struc¬ 
ture or other identifiable particles was called macroporosity. 
But in the context of numerical modeling those terms can 
take on different meanings, because there are three different 
ways to model void space. In the first way, which is appro¬ 
priate when the void space is very small compared to any 


length scale of interest, the porosity is considered a contin¬ 
uum material property and modeled as part of the equations 
of state. That continuum approach to porosity modeling is 
discussed below. In the second way, some codes allow a 
single numerical cell to contain both material and void, and 
the resulting behavior is determined by a mixture theory. 
In both of these first two approaches, the porosity scale is 
smaller than a calculation cell. In the third way, the void 
can be so large as to encompass an entire numerical zone, 
and zones without material can be scattered throughout a 
problem domain in some defined way. The differences be¬ 
tween the approaches will depend on the length scale of the 
modeled phenomena compared to the length scale of the 
porosity. 

The response of an asteroid to an impact is strongly af¬ 
fected by the presence of porosity. In the outgoing shock- 
wave, a porous material can undergo significant permanent 
compression and become hot, which creates a significant 
energy sink ( Asphaug et al., this volume; Davison et al., 
2010). That effect will be included in any numerical model 
where the size scale of the porosity is smaller then the width 
of the initial outgoing compression pulse. Since it is typical 
in a code to numerically smear a shock over several calcu¬ 
lation zones, even the third of the above porosity models 
can model significant crush up and energy loss, depending 
on the resolution of the calculation. In large scale collisions 
(say between bodies of a size of 100 km), a shock wave can 
lead to compression of porous bodies even if they contain 
large (~ kilometer size) voids, as long as these voids are 
smaller than the relevant scale (e.g. the impactor size). 

In addition to the effects at the shock, porosity and the 
material’s resulting crushability can also have a dramatic 
affect on the entire cratering process. Rather than an exca¬ 
vation processes, an entire crater can be formed by a down¬ 
ward flow crushing the material beneath the crater floor 
(Housen and Holsapple, in preparation, 2014). 

Recent experiments and the scaling theory for the regime 
of cratering dominated by target porosity will be discussed 
in section 3. The different approaches proposed to model 
the various effects of porosity described above will be pre¬ 
sented in section 4.4.5. 

2.2 Strength 

The outcome of an impact into an asteroid, whether a 
crater or a disruption, will ultimately be determined by 
gravity and some strength measure of the material of the 
object. There are many measures of strength for a geologi¬ 
cal material, and, over the last decade or so, that variety has 
been identified and is often included in our scaling theories 
and in numerical calculations. Strength measures can in¬ 
clude tensile strength, compressive strength, shear strength, 
crush strength and others; each governs the ability of a ma¬ 
terial to withstand a different kind of stress state. In the 
usual continuum approach, each of those strengths is char¬ 
acterized by a different portion of a single ’’strength enve¬ 
lope”: a boundary defined in stress space between elastic 
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and inelastic (permanent) deformation. An additional part 
of the modeling (flow rules) then describes the nature of 
the inelastic deformation from flow or fracture. In addition, 
the prior shock history can modify or ’’damage” the mate¬ 
rial, and that also must be accounted for. Holsapple (2009) 
presented a review of strength theories appropriate for ge¬ 
ological materials. The interested reader can refer to that 
reference, in addition to further detail below (section 4.4). 

2.3 Late-stage processes 

Since the times of the early lunar studies it is been ob¬ 
served that the large craters and basins have a substantially 
lower depth to diameter ratio than the smaller ones. For 
the large craters, the slopes of the outer walls are typically 
well below the angle of repose expected for soils and rocks. 
These were judged to be puzzling because the angle of re¬ 
pose typically determines the static equilibrium slope an¬ 
gles of soils and rocks. Melosh (1979) proposed that the 
effect was due to the dynamic weakening of rock in the 
latter stages of crater formation by the action of acous¬ 
tic vibrations. Since then, many other calculators include 
a ”late-time” period of crater formation using rheological 
models that suppose the presence of that acoustic fluidiza¬ 
tion. Those methods are common today. An alternative 
viewpoint was presented by Holsapple (2004a), but his ap¬ 
proach has yet to be fully developed. 

These approaches are presented below in section 4.4.7. 
An application of a dynamic weakening model in the case 
of large-scale collisions on Asteroid 4 Vesta (Jutzi et al., 
2013) is discussed in the chapter by Asphaug et al. (this 
volume). 

3. SCALING LAWS 

As mentioned in the introduction, the experiments we 
can make on Earth are not at the size scale, gravity levels, or 
impact velocities of interest to most of solar system impact 
events. For that reason, the results of experiments in the 
laboratory must be extrapolated, often over many decades, 
to predict the results of impacts into asteroids. How does a 
10 km asteroid behave compared to a 10 cm lab sample? 

The physical assumption forming the foundation of 
modern scaling theories is that of a general ’’point source”. 
Any impactor has three fundamental independent measures: 
a radius a, a velocity U, and a mass density 5. Equally well 
the three independent measures can be taken as the diame¬ 
ter, momentum, and kinetic energy or any other three inde¬ 
pendent combinations. In any case, they contain the three 
independent units of length, mass and time. Then, when 
that impactor collides at high velocity with an asteroid, it 
sets up a highly dynamic event affecting a region much 
larger than the impactor size, and over a timescale much 
longer than that of the initial deposition of energy into the 
surface. Physically, an appropriate assumption is that the 
deposition of momentum and energy is instantaneous into 
a region of vanishing dimensions compared to any length 


scale of interest such as the final crater. Of course, that 
assumption cannot be made for very low speed impacts or 
other cases where a final crater may be only slightly larger 
than the impactor. 

For a point source, the governing impactor measures 
cannot retain a length scale or a timescale. From that as¬ 
sumption it follows that the individual values for the size, 
velocity, and mass density do not affect the outcome. In¬ 
stead there can be at most one single combination of those 
3 variables that ’’measures” the impactor. That measure is 
then used in conjunction with those defining the material 
behavior of the target object to develop the scaling theory. 

The earliest point source solutions for impacts simply 
assumed that the correct measure was the kinetic energy of 
the impactor, that defined what is now called ’’energy scal¬ 
ing”. Dienes and Walsh (1970) noticed in code calculations 
of impacts into metals that the same crater was obtained for 
different impacts having the same all M where p = 0.58. 
They did not identify that result as the signature of a point 
source, nor did they recognize that the general form is quite 
universal, but with different exponents depending on mate¬ 
rial type. The Z-model of cratering by Maxwell (1977) was 
essentially another early point source model, but again it 
was not identified as such, and it was only applied to the ge¬ 
ometry of the cratering flow field. Holsapple and Schmidt 
(1982) developed crater scaling from the assumption of a 
general point source, measured by what they called a ’’cou¬ 
pling parameter”. They showed that it must always have the 
power law form C = aU p ‘5 l/ , but the governing exponents 
p and v cannot be immediately predicted, because they de¬ 
pend in a complicated way on the material of the target. 
However, it was proved that 1/3 < p < 2/3 (Holsapple 
and Schmidt, 1982). They also applied the same measure to 
a variety of outcomes of cratering and determined definite 
interrelations between the power laws for all outcomes of a 
given event. 

Over the years the implications of that theory have been 
applied to many different impact outcomes from both ex¬ 
periments and numerical simulations, including crater size, 
crater formation time, shock wave propagation time, catas¬ 
trophic disruptions, ejecta characteristics from the impact, 
momentum transfer to asteroids, and others. From that vari¬ 
ety of applications it is been well determined that the point 
source predictions work surprisingly well, and, for moder¬ 
ately porous materials p ~ 0.4, while for non-porous ma¬ 
terials p ~ 0.55, with only small variations found. And 
in many cases, the use of that assumption leads to sim¬ 
ple power-law scaling for many features of interest. Well- 
known examples are the power laws for crater dimensions 
when a crater size is determined by the surface gravity 
(’’gravity regime”) or when it is entirely determined by a 
strength (’’strength regime”), power laws for ejecta amounts 
and velocities, for stress decay, and others. 

The reader interested in the details and the earlier appli¬ 
cations might begin with the scaling review article in Hol¬ 
sapple (1993). A review and several applications including 
the catastrophic disruption cases were presented in the Hol- 
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Fig. 1.— Cratering on rocky asteroids can be dominated by sur¬ 
face spall phenomena when gravity and the crater are small. The 
lower shaded region of this plot shows the extent of that region, 
it would include all craters possible on bodies of km size, and all 
craters smaller then 200 meters on Eros. 

sapple et al. chapter in the 2002 Asteroids III volume. 

In addition to the now well-known strength and gravity 
regimes for impacts, two new regimes for cratering have 
been recently introduced. First, Holsapple and Housen 
(2013b) defined a “spall” cratering regime for small craters 
on rocky bodies, those are dominated by the tensile spall 
strength, and not the more common shear strength that de¬ 
termines excavation craters. Those have been well known 
for explosive craters smaller than a meter or so in rocky 
targets on Earth and in laboratory experiments in compe¬ 
tent rocks (e.g. Gault 1973), although that regime has usu¬ 
ally been ignored for planetary applications. But at the low 
gravity on an asteroid, that regime can include much larger 
craters. The extent of that regime as a function of surface 
gravity is depicted by the shaded region in Figure |T] from 
Holsapple and Housen (2013b). For a body such as Eros 
(16km) all craters smaller than about 1 km are predicted to 
be spall craters. They would be flatter and shallower than 
excavation craters, and would eject blocks, not ejecta with 
more uniform smaller particles. These analyses are new and 
require more thorough experimental and numerical investi¬ 
gation. 

Secondly, Housen and Sweet (2013) presented experi¬ 
ments and scaling of impacts into highly porous materials. 
Large porosity also adds a new regime, as illustrated in the 
schematic shown on a plot of scaled volume versus the grav¬ 
ity scaled 7t2 parameter (Figure [2]). 

As indicated above, the point source scaling theory has 
also been applied to catastrophic disruptions (e.g. since 
Housen and Holsapple, 1990). Recently, Leinhardt and 
Stewart (2012) extended this approach to define ’general’ 
scaling laws that included collisions between gravity domi¬ 
nated bodies of comparable sizes. While this approach was 
shown to work well for some specific regimes ( Leinhardt 
and Stewart 2012), its general applicability still remains to 
be validated. 

A great number of experimental studies in either the cra- 




Fig. 2.— The scaled size Hy depends upon the gravity scaled 
size II 2 in different ways for different materials, (a) For a cohe¬ 
sionless soil with low or moderate porosity such as dry sand, stan¬ 
dard ’’gravity scaling” applies at all size scales, and for increasing 
size or gravity there is a reduction in cratering efficiency, (b) For a 
small or moderately porous material with cohesion, the cohesion 
is the dominant contribution to the shear strength at small crater 
size scales, so that the cratering efficiency is constant. But for in¬ 
creasing size, there is a transition to the gravity-dominated regime, 
(c) A cohesionless material with high porosity is gravity scaled at 
small sizes, but at large sizes is created by compaction, indepen¬ 
dent of gravity, so Ily approaches a horizontal asymptote. ( d) A 
material with both cohesion and high porosity potentially shows 
all three of the cohesion, gravity, and compaction regimes. 

tering or the disruptive regime have been performed since 
the publication of Asteroids III. This chapter is focused on 
the modeling of collisions; therefore we will not review this 
great amount of experimental work. We just give a few 
references to the interested reader noting that these experi¬ 
ments greatly contributed to our understanding of the colli- 
sional process at small laboratory scales and always give a 
necessary point of reference to test numerical methods (see 
Section 4.6). In particular, various kinds of target materials 
have been considered such as targets made of sintered glass 
beads (e.g. Setoh et al. 2010, Hiraoka et al. 2011), pumice 
(Jutzi et al. 2009), porous gypsum, sometimes admixed 
with small pebbles or glass beads ( Okamoto and Arakawa 
2009, Leliwa-Kopystynski and Arakawa 2014, Yasui and 
Arakawa 2011), dense (soda lime or quartz) cores and 
porous (gypsum) mantles ( Okamoto and Arakawa 2008) 
and other porous materials (e.g. Housen and Sweet, 2013; 
Nakamura et al., 2014). Moreover meteorites (e.g., Flynn 
et al. 2008, Kimberley and Ramesh 2011) have been used as 
well as 1-meter diameter granite spheres to check the effect 
of size on the cratering outcome (Walker et al. 2013). In 
some cases, scaling parameters allowing an extrapolation at 
larger scales were also investigated. Holsapple and Housen 
(2012) study experiments and scaling for the momentum 
coupling of impacts as related to the deflection of Earth- 
threatening asteroids by impacts (see also section 5.2.). A 
comprehensive review of experiments and ejecta scaling has 
been published by Housen and Holsapple (2011). 
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The increasingly more sophisticated material modelling 
in shock physics code allows "numerical experiments” to 
be conducted for a much larger parameter space than it is 
possible to cover experimentally. Moreover the effect of in¬ 
dividual properties such as friction and porosity can be in¬ 
vestigated in detail. First results of this promising approach 
for the cratering regime have been presented in Wunnemann 
et al. (2010) and are summarized in section 5.1.1. Recent 
numerical studies of asteroid disruptions will be discussed 
in section 5.1.2. 

4. NUMERICAL MODELING 

4.1. Introduction 

As a complement to experimental and theoretical ap¬ 
proaches, numerical modeling has become an important 
component to study the outcome of collisions. In the last 
decade, advances in both computational resources and nu¬ 
merical methods have allowed the properties and processes 
involved (see section 2) to be modeled more and more real¬ 
istically. 

Here we present the principles of numerical impact mod¬ 
eling. We show some examples of codes used in the field 
and discuss the various model approaches. 

4.2. Numerical techniques 

The two most common approaches to simulate impacts 
and collisions use shock physics codes (’hydrocodes’) and 
particle codes. 

Shock physics computer programs use continuum the¬ 
ory, and can calculate the entire dynamical processes in¬ 
cluding the propagation of shock waves and resulting fields 
of displacement, velocity, strain, stress, etc., as function of 
time and position (e.g. Anderson, 1987). They rely on 
mathematical constitutive models: for the thermodynam¬ 
ics which often includes melt and vaporization, for the de¬ 
formation processes, and for the failure, fracture and flow. 
Those material equations are the outcomes of testing of ma¬ 
terials at the various states of interest. The equations are 
solved in a time-stepping manner on a geometrical com¬ 
putational grid (or interpolation points), usually with zones 
much smaller that the impactor dimension. The allowable 
time step size must be shorter than the time for the pas¬ 
sage of a shock wave across the smallest zone (’’Courant- 
Friedrichs-Levy” stability criterion). Depending on the nu¬ 
merical method used, additional time step restrictions are 
required (e.g. Anderson, 1987). It is not uncommon in 2D 
applications to include hundreds (n) of space zones in each 
direction, for a total of n 2 ~ 10 5 and runs for several 10 6 
time steps. In 3D problems, many more zones ( n 3 ) are re¬ 
quired, so generally n must be much smaller. That illus¬ 
trates the great advantage of codes that can calculate in 2D. 
However, many problems are inherently 3D. 

In contrast, particle codes (including ’’discrete element 
codes”) assume a collection of simple interacting particles. 


These codes generally only perform 3D calculations. Their 
collisions are modeled using heuristic descriptions for resti¬ 
tution, friction and viscosity, and their interactions include 
mutual inter-particle gravitational forces. The balance of 
linear and angular momentum is all that is required to cal¬ 
culate their motions and rotations, so there is no use of mass 
balance or energy balance concepts. They do not include 
crushing, melting, vaporization or other phenomena occur¬ 
ring in high-speed impacts, so these models are limited to 
low velocity, low stress events. The number of particles 
might be millions, but that is still many orders of mag¬ 
nitude less than the actual number in a soil-like structure 
in any problem of interest. Thus, it is inherently assumed 
that there are ’’enough” particles to approach the ’’infinite- 
number”, continuum limit. And, to date, it appears that the 
governing parameters need be chosen on a case-by-case ba¬ 
sis. 

There are two classes of particle models. The simpler 
and earlier approaches are the so-called ’’hard sphere” mod¬ 
els in which particle collisions, which are assumed to oc¬ 
cur instantaneously, are predicted in advance and are gov¬ 
erned entirely by coefficients of restitution chosen by the 
user. Such a model has been implemented in the N-body 
hard-sphere discrete element code pkdgrav (Richardson et 
al., 2000), which has been used for different applications 
in Planetary Science (see Michel et al. this volume for 
an application to asteroid family formation) and which has 
been adapted to enable dynamic modelling of granular ma¬ 
terials in the presence of a variety of boundary conditions 
(Richardson et al., 2011). Other examples of such codes 
are polyhedral rubble piles codes ( Korycansky and Asphaug 
2006; Movshovitz. et al. 2012). 

The hard sphere approach is equivalent to the early at¬ 
tempts to model fluid mechanics by a collection of very 
sparse interacting atoms or molecules. It works well for 
rarefied gas dynamics and it is well known that the aver¬ 
aging across many particles leads to the classical contin¬ 
uum equations of perfect fluids. However, dense systems 
involving multiple collisions and enduring contacts require 
another approach, such as the ’’soft sphere” discrete ele¬ 
ment method. In this approach, the particles are allowed to 
have finite interaction times governed by elastic concepts, 
and in principle even static cases of enduring contact can 
be included relying upon modest penetration between par¬ 
ticles. The soft sphere approach was recently included in 
the code pkdgrav (Schwartz et al., 2012); it is also used in 
other codes (e.g. Sanchez and Scheers, 2011). Moreover, 
these discrete approaches can include bonding forces be¬ 
tween particles, as a first analog of cohesive materials (see, 
e.g., Schwartz et al., 2013, Sanchez and Scheeres, 2014). 
Such capability is also now available in the commercial 
finite-element code, LS-Dyna. 

But this soft-sphere extension generally comes at a price. 
In effect, although the soft sphere approach has the advan¬ 
tage of not requiring collisions to be predicted in advance, it 
comes at the expense of much smaller integration timesteps 
than the hard sphere approach, which can limit the integra- 
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tion timescale. That is similar to the time step restrictions 
for the continuum codes. On the other hand, because it can 
be implemented into codes, like pkdgrav , that are fully and 
efficiently parallelized, it is currently possible to follow the 
evolution of millions of particles over a fairly large range 
of conditions. Murdoch el al. (this volume) give a review 
of discrete element methods and continuum approaches ap¬ 
plied to the dynamics of granular materials at the surface of 
asteroids. 

The use of particle codes to simulate rubble pile colli¬ 
sions at low speeds (below the sound speed of the materials) 
was discussed in Asteroids III (Richardson et al. 2002). In a 
recent study, Ballouz et al. (2014a, b) used the particle code 
pkdgrav to simulate low-velocity collisions between rotat¬ 
ing rubble pile in order to measure the effect of the initial 
rotation of colliding bodies on the outcome. 

In situations where the whole process of a large-scale 
asteroid collision is investigated, including both the initial 
shock passages, heating, fragmentation and the subsequent 
reaccumulation of fragments, a hybrid approach is often 
used (e.g., Michel el al., 2001, Michel et al., this volume). 
There the fragmentation is computed with a shock physics 
code and the gravitational reaccumulation with a particle 
code. This approach is limited by the numerical resolution 
that fixes the minimum size of tractable fragments (typi¬ 
cally down to ~ 10 meters for simulations involving krn- 
sized bodies). 

In the fragmentation phase of a hypervelocity asteroid 
collision, self-gravity can usually be neglectecQ because 1) 
the overburden pressure is small compared to the ampli¬ 
tude of the shock wave and 2) the fragmentation timescale 
is much smaller than the reaccumulation timescale. The 
fragmentation timescale is given by the time it takes for a 
shock wave to travel trough the whole target ry ~ Rt/c s 
where R t is the target radius and c s a wave speed. On 
the other hand, gravitational reaccumulation proceeds on 
a time scale of Td yn ~ (Gp) -1 / 2 ~ 2200 s (for a bulk 
density p = 3000 kg/m 3 and the gravitational constant G = 
6.67 x 10 -11 m 3 fc<? _1 s -2 ). For small bodies ( R t < a few 
100 km), Tf -C Tdyn and gravity does not affect the dy¬ 
namics of the fragments during the fragmentation phase. 
However, it is important to note that for very low veloc¬ 
ity collisions, gravity has to be computed during the whole 
process, even for small bodies. This is typically the case 
in accretionary collisions (see Asphaug et al., this volume) 
where the impact velocity v imp is of the order of the mu¬ 
tual escape velocity v esc = y/2G(M p + M t )/(R p + R t ), 
where M p is the mass of the projectile, M t the mass of the 
target and R p the projectile radius. For 10 (100) km diam¬ 
eter bodies, v esc ~ 5 (50) m/s (assuming a density of p = 
3000 kg/m 3 ). 

The hybrid hydrocode-particle code approach is detailed 
in Michel et al. (this volume), where also the newest colli¬ 
sion models in particle codes are presented. Here we focus 


1 However, for large asteroids it may be important to include fracture shield¬ 

ing due to compression 


on shock physics codes. 

4.3. Basic equations 


Shock physics codes solve the system of partial differ¬ 
ential equations that describe the conservation of mass, mo¬ 
mentum and energy for a continuous, compressible medium 
(see e.g., Collins et al., 2013 for a recent review). Examples 
and a description of such codes used in the field (SPH, CTH, 
iSALE, SOVA) will be given in section 4.5. 

The stress tensor is symmetric and is often divided into 
isotropic (hydrostatic) and deviatoric parts 

a ij = S ij - P6 ij (1) 


where the pressure P is the hydrostatic pressure, S IJ is the 
(traceless) deviatoric stress tensor and d , ' J is the Rroneker 
symbol. Using a Lagrangian reference frame, the conserva¬ 
tion equations (mass, momentum and internal specific en¬ 
ergy) can then be written using an indicial summation con¬ 
vention as follows: 



dp dv l 

dt +P d^=° 

(2) 


dv i 1 da ij i 

dt p dxi P 

(3) 

dE 

dt 

Pd 1 . .. 

= ——u* + -5‘ j r j 

p OX 1 p 

(4) 


where d/dt is the Lagrangian time derivative, p the density, 
v the velocity, E the specific internal energy, x the position 
and e is the deviatoric part of the strain rate tensor. The term 
g on the right hand side of eq. ([3]) accounts for any external 
acceleration, for instance due to gravity forces. 

To complete the set of equations, equations describing 
the material response are required. For the hydrostate com¬ 
ponents an equation of state (EOS) is required, which re¬ 
lates pressure, density and internal energy (or temperature, 
if its inclusion is convenient). And constitutive equations 
relating the deviatoric components are used. Both of these 
can be very complex, but have been developed over the 
years, often by the military community. The specifics are 
discussed in section 4.4. 

For many problems, the acceleration due to self-gravity 
is important and has to be taken into account in eq. ([3]). The 
components of the gravity acceleration can be computed by 
solving the Poisson equation for the gravitational potential. 
In particle based approaches (including SPH), the gravity 
acceleration for a particle k is directly given by 


g k = -Gj2 

k^l 


mi r kl 

r h r kl 


( 5 ) 


where G is the gravity constant, m the mass of the particles 
and r = r] the distance between the particles. For N par¬ 
ticles, the direct summation method leads to a complexity 
of 0(N 2 ) and is only practical when computers which are 
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designed to solve such problems (e.g. GPUs) are used. A 
method often applied for the self-gravity computation is the 
Barnes-Hut tree algorithm (Barnes and Hut, 1986) which 
allows to reduce the complexity to O (N log N). 

4.4. Material models 

An equation of state (EOS) relates density, internal en¬ 
ergy and pressure, and may also include porosity effects. 
A stiffness and strength model is needed to determine the 
deviatoric stress due to strains and possible material failure 
(section 4.4.2 and 4.4.3). 

Various processes occur during an impact on an asteroid 
and these have to modeled in a suitable way. Those may 
include the following: 

• effects due to porosity and other non-linearities 

- energy dissipation by compaction 

- damping of shock wave 

- reduction of wave speed in highly porous mate¬ 
rials 

• (post-impact) flow of granular material in the case of 
shattered asteroids or rubble piles 

• dynamical state weakening (e.g. important to model 
the collapse of large craters) 

4.4.1. Equations of state 
An often used form for the EOS is 


Lauson, 1972; Melosh 2007). In this model, the thermody¬ 
namics variables are derived from the Helmholtz free en¬ 
ergy. ANEOS includes a more accurate treatment of both 
melting and vaporization than the Tillotson approach and 
allows for other polymorphic and liquid/solid phase transi¬ 
tions. 

There are also complete tabular databases such as the 
SESAME library, developed by the Livermore National 
Laboratory. Those are efficient to use, valid over vast den¬ 
sity ranges, and commonly include complete thermody¬ 
namics including phase changes of melt and vapor. Often 
those are tabular compilations of the analytical forms. 

4.4.2. Strength models 

The strength model is fundamental for modeling impacts 
and collisions involving small bodies. It determines the ’im¬ 
pact strength’ in disruptive collisions, the size, final shape 
and characteristics of the crater in an impact, and so on. 

It is important to note that, depending on the loading 
conditions, various forms of ’strength’ exist (section 2). 
Here we give an overview of some strength models that are 
included in impact codes used in the field. 

A simple strength model commonly used is the von 
Mises model developed for ductile metals. In this model, 
plastic flow occurs at stresses larger than a single constant 
yield strength Yq. The von Mises criterion is implemented 
in shock codes by reducing the deviator stress (see. e.g. 
Benz and Asphaug 1994, 1995) by: 

S ij ■» fS ij (7) 

where / is computed by 


P = P{Pi E) (6) 

This form is convenient because in most hydrocodes, the 
specific internal energy (rather than the temperature) is 
computed directly. In fact, temperature is not required to 
solve the conservation equations. However, it may be use¬ 
ful to account for phase transitions and the thermal soften¬ 
ing when deviatoric stresses are included (section 4.4.2). 

The most simple equations of state have no thermody¬ 
namic coupling and the pressure is solely a function of den¬ 
sity (e.g. the Murnaghan EOS). A more sophisticated and 
widely used analytical equation is the Tillotson EOS, which 
was derived for high-speed impact computations ( Tillotson , 
1962). One of the major advantages of this EOS is its ef¬ 
ficiency. However, the Tillotson equation of state does not 
provide information about how to compute the temperature 
or the entropy of a material. Furthermore, the treatment 
of vaporization is not very sophisticated. For these reasons, 
the Tillotson EOS is mostly used to study impacts involving 
specific energies which do not lead to significant melting or 
vaporization. 

A more complex and thermodynamically consistent ana¬ 
lytical EOS model is ANEOS (M-ANEOS) (Thompson and 


f = min 


r y 2 
1 0 

.3 J 2 


,1 


( 8 ) 


This is commonly called the ’’radial return method” and is 
also, for the von Mises case, the direction of the ’’associated 
flow rule” of plasticity theories. Here, the second invariant 
of the deviatoric stress tensor J 2 = 1 ,S n -J S'4 is used as a 
scalar measure of the maximal shear stress. 

Although this model was developed for ductile materi¬ 
als, it was commonly used in impact calculations in geo¬ 
logical materials in the past. In combination with a tensile 
fracture model (see section 4.4.3) it gives reasonably ac¬ 
curate results in disruptive collisions. Jutz.i (2014) gives a 
comparison between this approach and more sophisticated 
strength models in simulations of disruptive collisions. 

However, common geological materials such as soils, 
rocks and ices have more complex behavior than ductile 
metals. An important characteristic of their strength is a 
substantial increase of shear strength with increasing con¬ 
fining pressure. That is the common feature of more so¬ 
phisticated failure criterion used for geologic materials. The 
general form for geological strength models is as indicated 
in the Fig. 1 of Holsapple (2009). There is a region in 
’’shear-pressure” space delimited by a closed curve at which 
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’’failure” or ’’flow” can occur. The ’’shear” can either be 
measured by the maximum shear stress r on any plane, or 
by the average shear given by the stress invariant The 
’’pressure” can as well be the actual average normal stress 
P or the maximum normal stress er on any plane. The 
Drucker-Prager model uses P and \/ d-i , while the Mohr- 
Coulomb model uses a and r. For small pressure, failure 
occurs at the upper limit pressure-dependent shear envelope 
defining the largest admissible shear. For large pressures, 
that shear envelope flattens. Then, at the right is a con¬ 
straint for the maximum pressure where compaction can oc¬ 
cur with compressive pressure. This limit is particularly im¬ 
portant for porous material for which ’’crushing” can start at 
very low pressures. Models for that right-most porous limit 
are described below. 

The initial slope of the shear envelope at low or tensile 
pressures is commonly called the ’’friction coefficient”, be¬ 
cause the form is similar to that used for friction between 
solid sliding blocks. But, in fact that is a misnomer, the 
actual mechanism is a result of the fact that a pressure im¬ 
pedes the movement of irregular-shaped grains up and over 
each other in a shearing flow. That slope then reduces at 
larger pressures. 

The Mohr-Coulomb and Drucker-Prager models include 
that pressure dependent yield strength and are the simplest 
failure models commonly used in soil and rock mechanics 
(e.g. Holsapple , 2009). 

One strength model (Collins et al., 2004) uses the fol¬ 
lowing pressure dependent shear strength for the intact ma¬ 
terial Yj : 


Yi = Y 0 + 


_ PiP _ 

1 + piPj(Y m - Y 0 ) ’ 


(9) 


where Yq is the shear strength at P = 0 and Ym is the shear 
strength at P = oo and p, L is related to the coefficient of in¬ 
ternal friction for the intact material. Here the ’’shear” is 
measured by \[J%. This model is implemented in a num¬ 
ber of codes (e.g. Collins et al., 2004; Senft and Stew¬ 
art, 2007; Jutzi, 2014). Other nonlinear functions includ¬ 
ing powers and exponentials are common (e.g. Hoek and 
Brown, 1980). 

To describe the yield strength of fully damaged rock, 
which includes granular material. 


Y d = p d P (10) 

is usecj^l where p d is related to the coefficient of friction of 
the damaged material. 

In addition, the model must specify how a stress exceed¬ 
ing the limit is mapped back to the failure surface in each 
time step (the ’’return method”) and also the resulting plas¬ 
tic flow increment during the time step (the ’’flow rule”). 
The yield strength Y is often used to simply reduce the de- 
viatoric stress (represented by sph) by a factor of Y/pp. 


2 For modelling regoliths in low gravity environments, it is not uncommon 
to add a (small) cohesion 


That is again the ’’radial return method”. The plastic strain 
increment is often assumed to be in shear only (no volume 
change). But in general plasticity theories, the ’’associated 
flow rule” is more common, it is assumed to be in the di¬ 
rection perpendicular to the failure envelope. That case 
includes the important phenomena of the dilantancy of a 
granular material when flowing. Such details are not yet 
included in most planetary code simulations, although such 
an approach has been implemented recently in the iSALE 
code ( Collins 2014). Much more study and use of such 
models is warranted. 

One might note the occurrence of two characteristic 
strength measures in such equations. The strength Y () at 
zero pressure is commonly called the cohesion. It may be 
zero for dry sands, or on the order of a few kPa to several 
MPa or more for cohesive materials, and even hundreds of 
MPa for small solid rocks. Then the shear strength increases 
with pressure to a maximum value of Ym, a characteristic 
of the strength of individual grains, which may be as much 
as a few GPa. The importance of these values depends on 
the other pressure scales in the problem, and especially on 
the gravity-induced lithostatic stress pgh at a depth h. If 
that gravitational stress is much larger than the cohesion, 
but still less than Ym, then the cohesion can be ignored, but 
a dependence on the angle of friction will still occur. That 
is the case in what is called the ’’gravity regime” of crater¬ 
ing which holds for instance for large craters on Earth. In 
the gravity regime, there remains a dependence on that an¬ 
gle of friction, but not on any other strength measure. That 
dependence will emerge from any calculation including the 
strength equation. But for extremely large events, the grav¬ 
ity stress will be larger even than Ym- For example, for 
collisions between planetary-sized bodies, that will be the 
case. Therefore, in giant impact simulations, for example 
to study the Moon forming impact (e.g. Benz et al., 1989; 
Canup and Asphaug, 2001; Reufer et al., 2012), any form 
of strength is typically ignored. 

There are additional factors which must be accounted for 
in the construction of any strength envelope. It is not static, 
but changes according to the present state of the material at 
any point. For example, it can depend on the instantaneous 
strain rate, or on the temperature, and can change size and 
shape as the material is strained. The community is just 
beginning to include all known effects in our mathematical 
models. 

For example, to account for the thermal softening of ge¬ 
ologic materials for extreme impacts, Y may be further re¬ 
duced according to: 

y^ytanhj^ 7 ^*-i)| (11) 

(Collins et al., 2004), or using an similar function of specific 
internal energy E: 

y-y(iG2) 

\ &melt J 

where E is the specific internal energy and E me i t the spe- 
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cific melting energy. It is important to note that the temper¬ 
ature and energy at which melt occurs are strongly pressure 
dependent and not simply a constant, and that fact should be 
included in any thermal softening equation. A possible ap¬ 
proach is to use the location of the present thermodynamic 
state compared to the melt boundary in the EOS space to 
soften the material. 

4.4.3. Fracture (or damage) models 

There is an additional complexity when a rock fractures. 
When a rock is stressed to failure, it breaks and becomes 
granular; then the failure envelope should ultimately change 
from one with cohesion to one appropriate for a fully gran¬ 
ular material such as a dry sand or gravel. Some approaches 
in the past did not do that, but simply assumed zero strength 
(as water) when damaged. A scalar damage parameter D 
describing the accumulation of tensile and/or shear frac¬ 
tures and/or pore crushing from undamaged (D = 0) to 
totally damaged (D = 1) is often used to interpolate be¬ 
tween the intact (eq. |9]» and the damaged (eq. p~0] > strength 

Y = (D-l)Y i + DY d (13) 

where Y is limited such that Y < 1*. The damage param¬ 
eter is computed using a fracture model and it may also be 
related to porosity models (e.g. Jutzi et al. 2008). 

Continuum fracture models (e.g. Grady and Kipp 
(1980)) often use an underlying structure with a preex¬ 
isting Weibull distribution ( Weibull , 1939) of cracks that 
grow and coalesce under tensile loading. In these models, 
cracks grow at a fixed speed c g once a flaw becomes ac¬ 
tive. This finite speed of crack grow naturally leads to a 
rate dependent failure of the material, as it is observed for 
rocks (see, for example, Housen and Holsapple, 1990 for 
the application to disruptions). In Asteroids III ( Asphaug et 
al. 2002), various aspects of the rate and size dependence 
fracture models were discussed. The implementation of the 
fracture models in shock physics codes are discussed in e.g. 
Benz and Asphaug (1994, 1995), Collins et al. (2004) and 
Jutzi et al. (2008). These models are generally overlaid in 
addition to the strength envelopes discussed above. 

Damage can also increase due to pore crushing during 
the compression of a porous material. This effect can be 
taken into account using a (linear) relation between disten¬ 
tion a and damage D (e.g. Jutzi et al. 2008). 

Damage is usually treated as a scalar parameter. How¬ 
ever, others have introduced tensor measures of damage 
(e.g. Lubarda and Krajcinovic, 1993). Due to their com¬ 
plexity those have not yet made their way into impact stud¬ 
ies, although tensor damage has been included in some re¬ 
cent codes (e.g. Owen , 2010). 

4.4.5. Porosity models 

In section 2 we discussed the various scales of porosity 
and the effects of porosity during an impact, such as the 


absorption of energy due to compaction, or the reduction of 
the wave speed. 

Depending on the scale of porosity, it can be modeled 
as either macroscopic voids (macroporosity) or by using a 
continuum, sub-resolution material model (microporosity), 
or as a combination of both as discussed in section 2. Here 
the size scale for demarcation is a computational cell size. 

A widely used continuum model is the ”P-alpha” model 
(Herrmann, 1969; Carroll and Holt, 1972), which uses the 
distention a defined by 

a = p s /p (14) 

where p is the bulk density of the porous material and p s 
is the density of the corresponding solid (matrix) material. 
A crucial assumption in this model is that it is the density 
and specific internal energy of the matrix material that de¬ 
termines the pressure. That is true when ignoring any en¬ 
ergy content related to the porosity such as surface energy. 
Therefore, the bulk pressure of the porous material P is then 
related to the pressure in the solid component (matrix) P s : 

P = -P s (p s , E s ) = -P s (ap, E). (15) 

a a 

A significant feature of this model is that any existing EOS 
for a nonporous material P s (p s ,E s ) can be used as the solid 
component of a porous material of the same composition. 

A compaction model is required to determine the his¬ 
tory of distention as a function of the history of pressure 
change. It is often assumed that the crushup is independent 
of shear stress, although Jutzi et al. (2008) relate the chang¬ 
ing rate of the distention and of the deviatoric stress tensor. 
The form of the model for the historic P-a Herrmann form 
includes a crush curve a = a c (P), at which pressure in¬ 
creases will decrease the distention, defining the crush-up of 
the material. In addition, the porous model defines unload¬ 
ing elastic curves. Those apply until the stress state again 
reaches a failure boundary. Herrmann (1969) does not dis¬ 
cuss that feature. The curves for loading and unloading are 
obtained by pressure tests that load and unload the mate¬ 
rial and are defined in the model with appropriate algebraic 
forms. The quadratic form assumed by Herrmann is not 
well suited to geological materials, but is easily changed. 
For instance, Jutzi et al. (2008, 2009) use a combination 
of two power-law functions to successfully reproduce the 
crush curve of pumice (Fig. 1 in Jutzi et al., 2009). 

Since the EOS has the underlying density as a function 
of both the pressure and the internal energy, the crush curve 
should include that also, so Holsapple (2008) suggests that 
a = a c (P, E). In particular, the crush strength should be¬ 
come zero when the thermodynamic state is at melt. That 
feature is not yet common in models. 

Elastic unloading only occurs until the stress state 
reaches some other point on the enclosing failure envelope. 
When shear stress is present, that can occur at positive pres¬ 
sure, but it will be at pure tensile pressure if the shear is zero 
and cohesion is present. At those states, increasing porosity 
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will occur: the occurrence with shear is the physical pro¬ 
cess of dilation in shear flow. Such an approach has been 
implemented recently in the iSALE code ( Collins 2014). 

The Herrmann approach to solving the history of disten¬ 
tion following these paths used a time sub-cycling which 
was essentially a forward differencing. The difficulties 
in that numerical method motivated an alternative model 
which determines the distention as a function of volumet¬ 
ric strain, it was presented by Wiinnemann et al. (2006). 
The so-called ’’epsilon-alpha” model addresses the above 
mentioned problem of the iterations in the ”P-alpha” model 
due to the interdependency of pressure P and distention 
a. When implementing eq. [15] and the crushing curve 
a = a(P) into a hydrocode for a given time step t, a t+ 1 
must be known to derive P t+ i, but a t +-\ = f(Pt.+ 1 ). A 
common solution is to use small sub-cycles to iterate the 
new P(t + 1) value. This method requires extra compu¬ 
tation time and may be numerically unstable under cer¬ 
tain circumstances. The problem is solved in the ’’epsilon- 
alpha”-model using the volumetric strain ey to determine 
the crushing of pore space: a = a(ey). The ’’epsilon- 
alpha” model distinguishes four compaction regimes, where 
the rate of compaction da/dey is calculated according to 
elastic compaction (ey < e e ), exponential compaction 
(e e < ey < e x ), power-law compaction (e x < ey < e c ), 
and the fully consolidated state (ey > e c ). 

The improved model ( Collins et al ., 2011) removed two 
shortcomings of the initial model. First, it accounts for the 
fact that the speed of sound of the pristine porous material 
can be substantially lower than the elastic wave speed in the 
solid material. A simple linear relationship for c(a) inter¬ 
polating between the speed of sound of the fully compacted 
material and the initial porous material is assumed. Sec¬ 
ondly, the improved version distinguishes between thermal 
and mechanical strains. The improvement uses an approx¬ 
imation about the equation of state, so further calls to the 
EOS subroutine are not needed. In contrast to the original 
’’epsilon-alpha” model the improved version is also appli¬ 
cable for highly porous material where material is heated 
extremely due to the compaction of pore space (PdV-work) 
resulting in thermal expansion of the solid component. 

An alternative idea that retains the P-a form was out¬ 
lined in Holsapple (2008). It is based on a Newton-Raphson 
approach rather than the forward differencing and elimi¬ 
nates the numerical problems in the original formulation. 

4.4.6. Fluidization models 

In the opinion of some researchers, the constitutive mod¬ 
els describing the strengths of rocks cannot explain the 
formation of complex crater structures with central peaks, 
where originally deep-seated material was uplifted several 
kilometers. And they do not explain the shallow outer wall 
slopes that are well below the angle of repose of the mod¬ 
els. It is assumed that an almost fluid-like rheology of mat¬ 
ter during crater formation is required, and it is from some 
weakening mechanism that lasts only temporarily. The lat¬ 


ter is an important constraint as steep, almost vertically 
standing flanks of central peaks can only be explained if 
rocks almost return to their initial strength. 

A successful approach to solve the problem of a tem¬ 
porary fluid-like rheology of rocks during crater formation 
was suggested by Melosh (1979). He proposed that heavily 
fractured and brecciated rocks behave like a granular flow 
excited by acoustic vibrations that have been interpreted to 
be observed in the ground motions generated by nuclear ex¬ 
plosions. The fluid-like nonlinear rheology decays as the 
amplitude of the vibration attenuates. The model requires 
that the wave-length of the acoustic signal is comparable 
to the size of the fragments. The original acoustic fluidiza¬ 
tion model has been simplified in the so-called block model 
where the acoustically fluidized material behavior is de¬ 
scribed by a Bingham model ( Melosh and Ivanov , 1999). 
The Bingham viscosity is proportional to the block size and 
the Bingham cohesion depends on the amplitude of the vi¬ 
bration. The latter is a function of time as the amplitude 
of the acoustic wave decays. A shortcoming of the acoustic 
fluidization or block model is the fact that neither block size 
nor decay time are know and can only be estimated for the 
size of a given structure. Wiinnemann and Ivanov (2003) 
suggested a heuristic linear relationship for both block size 
and decay time with the size of the impactor. They as¬ 
sume that larger impactors produce larger blocks and longer 
waves that attenuate slower than in case of smaller impacts 
where material is fractured into smaller fragments and the 
acoustic signal vanishes much quicker. The linear scaling 
of the block model parameters has been calibrated against 
the observed depth diameter ratios of the crater record of 
the moon (Wiinnemann and Ivanov, 2003) and on satellites 
(Bray et al., 2014). 

Figure [3] shows a simulation of the gradual transition 
with increasing crater size from simple to complex crater 
morphology, using that model. The parameter S is given 
by the ratio between strength Y and the hydrostatic pres¬ 
sure at maximum crater depth d max ■ S = Y/(pgd max ). 
The diagram illustrates the increasing crater depth while the 
transient cavity is growing and the decrease of crater depth 
when the crater floor collapses and starts to rise forming 
a complex crater with a central peak. According to some 
simple estimates by Melosh (1979) the critical value for 
crater collapse is given by S = 0.25 which corresponds 
to strength of the rocks of a few MPa, approx, one order 
of magnitude smaller than typical values for rocks. The 
application of the block model allows to assume realistic 
strength values for crater collapse at the observed transition 
diameter between simple and complex crater morphology 
on different planets. 

The approach described above is well established and 
was successfully applied in the past in many studies to 
model the formation of complex craters. However, there 
is an ongoing debate about the actual underlying physical 
mechanism that causes the temporary weakening. Holsap¬ 
ple (2004a, 2013) suggested that existing strength theories, 
in their complete form, can model the late time readjust- 
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Fig. 3.— Crater depth versus time for different impact condi¬ 
tions (S = 1.02-4.46). The length parameter is normalized by the 
maximum crater depth dmax, time is scaled by the ratio of impact 
velocity v t and dmax- The no-fluidized rheology is given by a 
Drucker-Prager yield strength envelope where Y = Yq + <j>p, with 
To=25 GPa and <j>= 0.1. Examples of the resulting crater morphol¬ 
ogy for different S-values (different projectile sizes) are shown 
above and below the diagram. Note, all length scales are normal¬ 
ized by d m ax (taken from Wiinnemann and Ivanov , 2003). 


ments discussed above, but his approach has not yet been 
fully vetted. Some important evidence for the related prob¬ 
lem of the collapse and runout of a granular cliff is given in 
section 4.6.3 below. 

4.5 Examples of hydrocodes 

Various approaches are used in shock physics codes to 
solve the continuum system of partial differential equations 
detailed in section 4.3. The equations may be cast in a “La- 
grangian” reference frame that follows the material or in 
an “Eulerian” reference frame that is fixed in space. The 
continuum codes include finite difference, finite element, 
and SPH methods. Each numerical method has its own 
strengths and weaknesses. In a recent review by Collins 
el al. (2013), the various methods are presented and numer¬ 
ical issues like resolution, the treatment of shocks, multi¬ 
material approaches, etc. are discussed. Here we present 
a few examples of commonly used hydrocodes. For each 
method, the specific material models used and the pros and 
cons are indicated. 

4.5.1. Grid based codes 

Most numerical simulations of planetary collision pro¬ 
cesses use Eulerian grids, because large deformations are 
difficult to track in the Lagrangian approach. 

They typically use a two-step approach, where the first 
step is the Lagrangian step where the deformation of the 
grid according to a given velocity field is calculated, and 
then there is a second step that maps the grid back on its 


original location in space. The remapping requires special 
treatment of material boundaries which is usually solved by 
tracking or reconstructing the interface between different 
matter (for a discussion of different methods of interface 
tracking/reconstruction see Elbeshausen and Wiinnemann, 
2010). The main advantages of grid-based codes are: 

• Lagrangian, Eulerian or an arbitrary spatiotemporal 
transition between both reference frames is possible 

• High resolution that may be increased locally by the 
adaptive mesh refinement (AMR) method (Berger 
and Oliger, 1984; Berger and Colella, 1989) em¬ 
ployed in CTH 

• Adequate treatment of solid state deformation, liquid 
flow, and gas expansion including appropriate rheol¬ 
ogy models 

• Coupling of the grid-based methods with massless 
Lagrangian tracer particles to track the spatiotempo¬ 
ral thermodynamic history of matter 

• Relatively straight forward implementation of a vari¬ 
ety of Equation of State (EOS) models 

Grid-based models come also with some disadvantages: 

• In Eulerian mode, material interfaces are not tracked 
perfectly and bulk properties must be defined for 
mixed cells 

• Matter is treated as continuum which makes it diffi¬ 
cult to treat fragmentation properly. Separate frag¬ 
ments tend to be underresolved unless AMR is em¬ 
ployed 

• Although multi-material handling is possible, mostly 
the codes do not deal with multi-phase processes 
where each phase moves at its own speed resulting 
in mixing of matter. Such processes are in particular 
important for the interaction of solid state fragments 
with an expanding vapor plume. 

In this paper we mention three shock physics codes 
widely used in planetary science. A much more complete 
overview of available hydrocodes is given in Pierazzo et al., 
(2008). 

The most advanced software package among those is 
the Sandia National Laboratories CTH code (McGlaun et 
al., 1990); However, because of its military uses, this code 
is limited in its availability to US citizens. Besides highly 
advanced material models it also includes AMR and self¬ 
gravity. For disruptive large scale collision processes such 
as the moon-forming impact event self-gravity is essential 
('Can up et al., 2013). Canup et al. (2013) also compare re¬ 
sults produced by different approaches such as grid-based 
models (CTH) and mesh-free codes (SPH, see below). 
Other widely used hydrocodes are SOVA (Shuvalov, 1999) 
and iSALE (Amsden et al., 1980; Collins et al., 2004; 
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Wiinnemann et al, 2006; Elbeshausen et al., 2009) that 
are, with some limitations, freely available for scientific 
purpose. SOVA contains a specific routine to deal with the 
interaction of lithic and molten ejecta and the vapor plume. 
As mixing of all phases occurs it may be better called ejecta 
plume. This process requires so-called multi-phase hydro¬ 
dynamics where at least two different phases (ejecta and va¬ 
por) travel at different speeds interacting and mixing with 
each other. SOVA addresses this problem by introducing 
representative tracer particles where each tracer represents 
a certain number of fragments of a given size. The tracers 
exchange momentum and energy with the surrounding gas, 
but not with each other. The approach only works if the vol¬ 
ume of fragments is small relative to the volume of the gas. 
The method has been successfully applied to model dusty 
flows ( Shuvalov , 1999), tektite formation and deposition 
(,Stoffler et al, 2002), and the development of Chicxulub 
distal ejecta ( Artemieva and Morgan, 2009). 

The iSALE hydrocode is based on the SALE hydrocode 
solution algorithm ( Amsden et al., 1980). To simulate 
hypervelocity impact processes in solid materials SALE 
was modified to include an elasto-plastic constitutive 
model, fragmentation models, various EoS (Tillotson and 
ANEOS), and multiple materials ( Melosh et al, 1992; 
Ivanov et al., 1997). More recent improvements include a 
modified strength model ( Collins et al, 2004) and a poros¬ 
ity compaction model (Wiinnemann et al, 2006; Collins 
et al., 2011). The 3D version uses a numerical solver as 
described in Hirt et al. (1974). The development history 
of iSALE-3D is described in Elbeshausen et al. (2009). 
iSALE has been used to model the collision of highly 
porous planetesimals (Davison et al, 2012), impacts on 
the surface of Vesta ( Krolin et al, 2014) and the formation 
of large craters on Lutetia (Cremonese et al, 2012), and the 
effect of an oblique impact angle on crater formation (Elbe¬ 
shausen et al, 2013; Collins et al, 2012; Elbeshausen et 
al, 2009). 

4.5.2. Smoothed Particle Hydrodynamics (SPH) codes 

Smoothed particle hydrodynamics (SPH) is a meshless 
continuum Lagrangian approach (and in spite of its name is 
not a ’’particle” code). In the basic SPH method, approx¬ 
imate numerical solutions of the fluid dynamics equations 
are obtained by replacing the fluid with a set of particles 
(see e.g., Gingold and Monaghan 1977, Monaghan 2012). 
The properties of these particles are smoothed over a cer¬ 
tain length h by a kernel function. That provides the link to 
a continuum approach. As the flow evolves, material mass 
moves with the particle and local density is calculated based 
on the proximity of nearby particles. The main advantages 
of this mesh-free method are: 

• SPH is a very robust scheme. It is ’straightforward’ 
to add new physics. 

• Large deformations of the simulated objects are han¬ 
dled easily. 


• The description of free surfaces is trivial. 

• Due to the particle nature of the method, fragmenta¬ 
tion is modeled naturally by the separation of particle 
clusters. 


Disadvantages of the SPH method include: 

• The spatial resolution is typically lower than in grid- 
based simulations. 


• Boundary conditions and discontinuities (e.g. mate¬ 
rial interfaces) are complicated to handle. 

• It is mostly applied in 3D, which often requires com¬ 
putationally expensive million-particle calculations 


In the past years, the range of applications of the SPH 
algorithm has increased significantly. A number of SPH 
codes have been recently developed to study problems in 
planetary sciences. Examples are a SPH code for the 
modeling of pre-planetesimals (Geretshauser et al. 2011), 
SPHERAL (Shapiro, et al. 1996, Owen 2010), GADGET 
(Springel, 2005; Marcus et al 2009), SPHLATCH (Reufer 
et al., 2012) or a GPU-based SPH code (Kaplinger et al. 
2013). 

A widely used SPH code to model collisions among 
rocky bodies was developed by Benz and Asphaug (1994, 
1995). This code was further extended by Jutzi et al. (2008, 
2013) with the goal to realistically model rocky bodies with 
various internal structures (see section 4.5 for a comparison 
to laboratory experiments). The most recent version of this 
code (see Jutzi, 2014) includes a pressure depended strength 
model as outlined in section 4.4.2, a tensile fracture model 
(section 4.4.3) and a porosity model based on the P-alpha 
model (section 4.4.5) and self-gravity. Friction is modeled 
either using the Coulomb dry friction law (eq. 10 1 or the 
rate-dependent model suggested by Jop et al. (2006). Self¬ 
gravity (section 4.3) is included as well. Recent applica¬ 
tions of this code in asteroid studies are presented in section 
5 and in the chapter by Asphaug et al. (this volume). 


4.6 Comparison to laboratory experiments 


An important step in the development of numerical 
methods (which includes the implementation of complex 
material models) is the validation against laboratory exper¬ 
iments. As rigorous testing is limited by the availability of 
appropriate experiments, including measurements during 
the highly dynamic processes, a preliminary step is bench¬ 
marking of different codes against each other. An example 
for a very successful benchmark and validation study is 
given in Pierazzo et al. (2008) where the modelling re¬ 
sults of specific test problems from 8 different hydrocodes 
were compared against each other and validated against 
laboratory cratering experiments, but only for water and 
aluminium. 

Here we provide some additional examples dealing with 
the fracturing and fragmentation resulting from collisional 
processes, and a granular flow problem. 
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Fig. 4.— Comparison between laboratory impact experiments 
(left) and SPH code calculations (right). The targets are pumice 
cubes of size of 7 cm with a porosity of ~ 70%. The impact ex¬ 
periments were performed at the Institute of Space and Astronau- 
tical Science (ISAS) of the Japan Aerospace Exploration Agency 
(JAXA) using a two-stage light-gas gun. Two impacts are shown 
at different times: t = 8 ms (top) and t = 1.5 ms (bottom). From 
Jutzi et al. (2009) 

4.6.1 Modeling the fragmentation of porous pumice 

Figure [4] shows a comparison of an SPH code calcula¬ 
tion (using the SPH method as described in section 4.5.2) 
to laboratory impact experiments. The target used in the 
experiments is porous pumice with a porosity of ~ 70%. 
The fragment size distributions resulting from the impacts 
with velocities ranging from 2 to 4 km/s could be well re¬ 
produced in the SPH code calculations (see Jutzi et al. 2009 
for details). This example illustrates the capabilities of the 
SPH method to model impact fragmentation. 

4.6.2 Fracturing in cratering experiments 

Fracturing caused by hypervelocity impact does not nec¬ 
essarily result in a complete disruption of the target. Fig¬ 
ure 0 shows the cross-section of the crater that was formed 
by the 5.4 km/s impact of a 1cm diameter iron projectile 
in ~20% porous sandstone. The experiment was carried 
out in the framework of the so-called MEMIN (Multidisci¬ 
plinary Experimental and Modeling Impact Research Net¬ 
work) project aiming at the validation of hydrocode mod¬ 
elling of hypervelocity impact processes (see Kenkmann et 
al., 2013 and papers therein). The left frame of Figure [5] 
shows a snapshot of the iSALE cratering model at 750 ps. 
The model includes the e — a porosity compaction model 
and a strength and damage model as described in (Collins 
et al., 2004). Some static strength parameters of the sand¬ 
stone were available (Kenkmann et al., 2011), others such 
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Fig. 5.— Comparison between laboratory cratering experi¬ 
ments (right) and iSALE simulation (left). The target was a 
100x100x50cm sandstone block with 20% porosity. The so-called 
Seeberger sandstone is composed of 97 w% Si02 with an uncon- 
hned compressive strength of ~ 60 MPa. The projectile was a 
lcm iron sphere impacting at 5.4 km/s (for further details see 
Kenkmann et al. (2011). Right: The thin gray lines indicate 
macroscopic flaws, the thick gray lines indicate crater profiles in 
different directions. Left: snapshot of the numerical model show¬ 
ing contours of the damage parameter. The black lines indicate 
the zone where tensile failure occurs and material spalls off the 
surface. Courtesy by the MEMIN-team 

as crushing strength of the porous sandstone were estimated 
or adjusted to match the observed crater depth. 

A more detailed calibration and validation of the crush¬ 
ing of pore space and sandstone is presented in (Giildemeister 
et al., 2013; Kowitz et al., 2013). The crater in the experi¬ 
ment is enlarged by spallation. Modelling of the actual spall 
of material from the surface is not included in the iSALE 
model; however, the spall zone is indicated by cells that 
have experienced failure. 

The damage zone in the model (red and orange contours 
represent the damage parameter as introduced above) un¬ 
derneath the crater shows some qualitative similarities with 
the flaws highlighted in red on the cross-section of the cra¬ 
tering experiment. On the scale of individual grains it was 
shown that the zone where pore space was crushed or frac¬ 
turing occurs in grains extends approximately to the same 
distance as in the model. 

4.6.3 Cliff collapse problem 

The cliff collapse problem is a useful test case for the 
pressure-dependent models used in code calculations, at 
least for low pressure. In that problem, an initial vertical 
’’cliff’ of material is suddenly released, and the mass falls 
and flows out to considerable lengths. 

This problem is of special interest for cratering because 
the same requirements of a temporary fluidization mecha¬ 
nisms have been asserted as necessary for the mechanics 
of landslides. The final angles of the runouts are typically 
< 10°, far below the angle of repose of geological materi¬ 
als. 

Holsapple (2013) presented detailed theoretical argu¬ 
ments and numerical calculations of that landslide problem, 
using only the standard Drucker-Prager soil model with a 
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35° angle of friction. He concludes that no additional ad- 
hoc models are required to reproduce observed laboratory 
results. Specifically, his numerical simulations correctly re¬ 
produce laboratory experiments with final slopes < 10°. A 
detailed investigation of those simulations explains the rea¬ 
son. It lies in the fact that the ’’slopes” allowed by the angle 
of repose argument are not violated, it is just that during the 
highly dynamic flows, the ’’slopes” must be measured rela¬ 
tive to the combination of the local gravity direction and the 
direction of the inertial forces. 

Jutzi (2014) also presented SPH calculations of that 
problem. He compared two models, one with a constant co¬ 
efficient of friction (eq. [T 0 | and one using a rate-dependent 
and a particle size dependent relation defined in an ’’inertial 
number” as suggested by Jop et al. (2006). It was found 
that both models reproduce the experiments very well. The 
two models lead to the same results because in this granu¬ 
lar flow regime, the inertial number stays small and there¬ 
fore the rate dependency is negligible. This finding provides 
further evidence that the global outcome of such events is 
well reproduced by using a simple Coulomb dry friction law 
with a single parameter pd, although at higher pressures a 
non-linear model is undoubtedly necessary. 

These examples show that the cliff collapse problem 
and the resulting runouts (at laboratory conditions) can be 
well reproduced with continuum codes using conventional 
rock mechanics. However, there is evidence that additional 
weakening mechanism (section 4.4.6) are necessary to ex¬ 
plain other cases, such as the very long runout slides ob¬ 
served on Mars (e.g. Lucas and Mangeney, 2007; Harrison 
and Grimm, 2003). 


5. APPLICATIONS 

In this section, some examples of recent modeling are 
presented which illustrate the various the effects of the ma¬ 
terial properties, target structures and impact conditions on 
the outcome of a collision. 

5.1. Scaling laws from numerical modeling 


5.1.1. Cratering regime 

Simple power-law relationships, so-called scaling laws, 
have been proposed to relate the size of the transient crater 
with the impact parameters such as the mass and velocity 
of the impactor (see section 3). Those have been shown to 
be a theoretical consequence of the point source assump¬ 
tion, as mentioned above. The scaling parameters in these 
equations (the so-called velocity exponent p and a propor¬ 
tionality factor that may be called I \) have been determined 
by laboratory experiments, mainly in sand or other granular 
materials. One way to determine the scaling parameters is 
to vary the so-called gravity-scaled size of an impact event 
7 T 2 = ( Lg)/v 2 , with the diameter L of the impactor, or the 
gravity g, or the impact velocity v, over a large as possible 




Fig. 6 .— Scaling parameters p and K 2 as a function of coeffi¬ 
cient of friction / and porosity <j>. (a) The scaling exponent p was 
determined by numerical impact experiments into material with 
different coefficient of friction, zero cohesion, and zero porosity. 
The dashed line was taken from Elbesliausen et al. (2009). (b) 
The scaling factor K 2 was determined by numerical impact ex¬ 
periments into material with different porosities, a coefficient of 
friction / = 0 . 8 , and zero cohesion. 
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Fig. 7.— Snapshots of the approximate transient crater for differ¬ 
ent impact angle and coefficient of friction. Relative change of the 
crater volume Vo normalized by crater volume for an equivalent 90 
degree impact V 90 with impact angle for (b) / = 0 (fluid-like tar¬ 
get rheology) and (c) / = 0.7 (sand-like target rheology). Note, 
the results from numerical models over a range of different 7 T 2 
values fall in between two lines calculated according to the verti¬ 
cal velocity component approximation (~ sin(4) 27 ; lower curve), 
and another line (~ sin($) 7 ). Taken from Elbeshausen et al., 
(2009). 


range. For example, one can increase the gravity using a 
geotechnic centrifuge, and decrease the gravity using drop 
towers. The geotechnic centrifuge method was developed 
both for explosive cratering and for hypervelocity impact 
cratering by Schmidt (1977) and by Schmidt and Holsapple 
(1980). Since gravity can be varied over two orders of mag¬ 
nitude, a 500 G centrifuge test models 1G events with a 500 
times larger scale. A meter-sized test at 500G has the same 
physics as a 0.5 km event at 1 G. However, there is no af¬ 
fect of the increased gravity at those magnitudes except for 
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materials with very little cohesion, so studies are primarily 
limited to granular materials such as sand. Although the 
properties of the target material such as porosity or inter¬ 
nal friction are somewhat under the experimenter’s control, 
varying these parameters independently is difficult. 

An alternative approach is to use numerical simulations 
to carry out numerical experiments which are not limited by 
experimental constraints in the laboratory. Of course, due to 
the numerous modeling questions mentioned above, those 
are not without their own uncertainties. And, as long as the 
point-source approximation is valid, they cannot contradict 
the theoretical scaling forms, but they can identify the un¬ 
known scaling coefficients. And they allow extensions of 
the parameter space from laboratory scale to natural dimen¬ 
sions under any given gravity regimes of planets, moons, 
and asteroids. The effect of material properties of the tar¬ 
get such as porosity, friction, cohesion can be investigated 
independently (Wtinnemann et al., 2010) and also the effect 
of the angle of incidence can be investigated ( Elbeshausen 
et al., 2009). In these studies it was shown that the velocity 
exponent p depends on the coefficient of friction / of the 
target material (Figure [ 6 ]» and the scaling factor K may be 
expressed as a function of porosity <j>. 

To incorporate the angle of impact 9 into the scaling 
equations it was suggested by Chapman and McKinnon 
(1986) and others based on impact experiments in sand to 
replace the impact velocity ve by the vertical component 
of the impact velocity r;j . In this case, for instance crater 
volume scales with the sine of the impact angle raised to the 
power 27 where 7 = 3p/(2+p). Extensive parameter stud¬ 
ies using models of impacts at different angles into granular 
targets with varying coefficient of friction (Figure [7j0 show 
that only for a typical coefficient of friction of / = 0.7 
this approximation holds true (Figure [ 7 }:, Elbeshausen et 
al., 2009). The vertical velocity component approximation 
also provides good estimates of crater volume and diame¬ 
ter for strength dominated craters in cohesive ductile targets 
such as metals (Davison et ah, 2011). For impacts in gran¬ 
ular targets with smaller coefficient of friction than typical 
for sand the size of the resulting crater is underestimated by 
this simple assumption (Figure [ 7 ] 5 ). Details of the impact 
may not be so easily matched using this assumption. 

5.1.2. Disruption regime 

To characterize the outcome of a disruptive collision, the 
critical specific impact energy Q* D which results in the es¬ 
cape of half of the target’s mass in a collision is often used. 
The parameter Q* D is called the catastrophic impact energy 
threshold (also called the dispersion threshold). The spe¬ 
cific impact energy is often defined as Q = 0.bm p v p /Mx, 
where m p , v p and Mr are the mass and speed of the pro¬ 
jectile and the mass of the target, respectively. The catas¬ 
trophic disruption threshold Q* D is then given by the spe¬ 
cific impact energy leading to a largest (reaccumulated) 
fragment M[ r containing 50% of the original target’s mass. 
In recent studies (e.g. Stewart and Leinhardt, 2009; Lein- 


liardt and Stewart 2012), a more general definition of the 
specific impact energy was proposed which also takes the 
mass of the impactor into account, which can be substan¬ 
tial in very low velocity impacts of near-equal-sized bodies. 
The corresponding radius Rci is then defined as the spher¬ 
ical radius of the combined projectile and target masses at a 
density of 1 g cm -3 . According to this new definition, the 
catastrophic disruption threshold is then called Q* RD - 

Values of Q* D (or Q* RD ) have been estimated using both 
laboratory and numerical hydrocode experiments (see e.g. 
Holsapple et al., 2002; Asphaug et al. 2002). For high ve¬ 
locity asteroid collisions, the first suite of numerical calcu¬ 
lations aimed at characterizing the catastrophic disruption 
threshold in both the strength regime and the gravity regime 
was performed by Benz and Asphaug (1999), who used a 
smoothed particle hydrodynamics (SPH) code Benz and As¬ 
phaug (1994,1995) to simulate the breakup of basalt and 
icy bodies from centimeters-scale to hundreds kilometers 
in diameter. More recently, Leinhardt and Stewart (2009) 
computed Q* D curves using the hydro code CTH McGlaun, 
1990 to compute the fragmentation phase and the N-body 
code pkdgrav to compute the subsequent gravitational evo¬ 
lution of the fragments. In this study, the dependency of 
Q* d on the strength of the target was investigated. In a re¬ 
cent study by Jutzi et al. (2010), the effect of target porosity 
on Q* d was investigated using an extended version of the 
SPH code (Jutzi et al., 2008). In this study, the size and ve¬ 
locity distribution of the fragments was computed as well, 
using the pkdgrav code. Benavidez et al. (2012) performed 
a study of a large number of collisions among R t = 50 km 
rubble pile bodies using the original SPH code by Benz and 
Asphaug (1994,1995). As discussed in section 3, Leinhardt 
and Stewart (2012) proposed general scaling laws for colli¬ 
sions among gravity dominated bodies. 

In the chapter by Asphaug et al. (this volume), a recent 
systematic study of the relative effects of various asteroid 
properties on the disruption threshold Q* D (Jutzi, 2014) is 
presented. 

5.2. Momentum transfer in small impacts 

The study of the momentum transferred in an specific 
impact on an asteroid as a function of impact conditions 
and the internal structure is crucial for performance assess¬ 
ment of the kinetic impactor concept of deflecting a poten¬ 
tially hazardous asteroid from its trajectory (see Harris et 
al. this volume). The momentum transfer is characterized 
by the so-called momentum multiplication factor j3, which 
has been introduced to define the momentum imparted to an 
asteroid in terms of the momentum of the impactor: 

P = 1 + Veil ( M p v p ), (16) 

where p e j is the escaping ejecta momentum, and M p and v p 
are the mass and velocity of the impactor, respectively (see 
e.g. Holsapple, 2004b). In the limiting case of an impact 
which produces no escaping ejecta, /? = 1 , and the mo¬ 
mentum transferred corresponds to the momentum of the 
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projectile (inelastic collision). However, in the case of an 
impact which produces a lot of material (ejected in the op¬ 
posite direction) with velocities larger than the escape ve¬ 
locity, we can have > 1 due to the contribution of p e j. 
Using the /3 factor, the resulting momentum of the target 
A P t (along the impact direction) can be obtained by 

A Pt = PPp, (17) 

where P p = rn p v p is the momentum of the projectile. 

It is important to note that only fragments with ejection 
velocities larger than the escape velocity of the target as¬ 
teroid, v e j ec t > v esc , escape the body and contribute to 
the momentum transfer. Moreover, the ejected material is 
slowed down during its escape and its trajectory is changed 
due to the gravitational attraction of the target. Therefore, to 
compute the escaping ejecta momentum, the value and an¬ 
gle of the velocity at infinity have to be used (see e.g. Hol- 
sapple and Housen, 2012; Cheng, 2012; Jutzi and Michel, 
2014). 

Holsapple and Housen (2012) presented scaling laws to 
extrapolate measurements of the momentum transfer in lab¬ 
oratory experiments to asteroid scales. 

Recent code calculations of the momentum transfer 
efficiency in impacts on asteroids for conditions typical 
for a kinetic impactor were performed by Holsapple and 
Housen (2013a) and Jutzi and Michel (2014). Holsapple 
and Housen (2013a) presented numerical simulations both 
for porous and non-porous materials, and obtained very 
good agreement with their laboratory results. 

In the study by Jutzi and Michel (2014), the effect of 
different degrees and scales of porosity on the momentum 
multiplication factor j3 was investigated for a range of im¬ 
pact conditions. 

Two kinds of porous target structures were considered 
by Jutzi and Michel (2014): (a) homogeneous: microp- 
orous only, and (b) heterogeneous: both micro- and macro- 
porous. For both structures, the microporous part of the 
material has a porosity of 50%. For structure (b), in ad¬ 
dition to the microporosity, macroscopic cracks with a size 
scale of ~ 0.3 m are randomly distributed in the target. The 
resulting total macroscopic void fraction is 10 %. 

The impact simulations were performed using a Smooth 
Particle Hydrodynamics (SPH) impact code (e.g. Benz and 
Asphaug 1994,1995; Jutzi et al. 2008, Jutzi, 2014) as pre¬ 
sented in section 4.5.2. 

In Figure [8] the outcome (in terms of damage) of an im¬ 
pact at 10 km/s is shown for the two target structures. 

In Figure[9] the results of the calculations of 3 are shown 
for the considered range of impact velocities (0.5-15 km/s) 
and the two target structures. At low impact velocities, the 
amount of momentum transferred is smaller using structure 
(b) than structure (a). However, these differences disappear 
at high velocities. This means that for porous targets, in¬ 
homogeneities at the scales considered here have negligible 
effects on the amount of transferred momentum for veloci¬ 
ties > 5 km/s. That confirms the view expressed earlier that 
the exact form of the modeling of porosity will not matter 



Fig. 8.— Left: damage (dark gray zones) produced by the impact 
at 10 km/s on a microporous target; the top figure shows the target 
from above, while the bottom figure shows a vertical slice. Right: 
same for a target containing both microporosity and macroporos¬ 
ity. From Jutzi and Michel (2014). 



Fig. 9.— Momentum multiplication factor /3 — 1 as a function 
of impact velocity for the two considered structures (a: homo¬ 
geneous microporous, and b: heterogeneous, micro and macrop- 
orous). From Jutzi and Michel (2014). 

for any phenomena with scales large compared to the size 
scale of that porosity. 

While the effect of target inhomogeneities on the mo¬ 
mentum multiplication factor /3 appears to be quite small, 
the effect of various material properties such as the tensile 
or crushing strength was found to be more significant (Jutzi 
and Michel, 2014). 

The results of this study verify that the momentum mul¬ 
tiplication factor j3 is small even for very high impact veloc¬ 
ities (/? < 2 for Vi mp < 15 km/s) in the case of porous tar¬ 
gets (with ~ 50% porosity). This is consistent with scaling 
laws (section 3) which, in combination with the results of 
laboratory experiments, predict a small value of /3 ~ 1 — 2 
for porous materials and comparable impact velocities (see 
Table 3 in Holsapple and Housen, 2012). It is also consis¬ 
tent with the numerical simulations performed by Holsap¬ 
ple and Housen (2013a). 

However, that is not true for non-porous targets. Val¬ 
ues of 3 as large as 5 have been directly measured at im¬ 
pact velocities of v,j mp < 5 km/s ( Holsapple and Housen, 
2013a) for rocky targets. Furthermore, the scaling predicts a 
marked increase like U 0 ' 5 for those targets, so at Vi mp < 20 
km/s those values might reach as much as /? ~ 10. If so. 
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that would significantly facilitate deflections by impacts. 

An important question is how to scale strength (tensile 
strength and crush-curve parameters) of a porous material 
to larger sizes. The strength properties of real asteroid ma¬ 
terials and their size dependency are not well constrained 
(e.g., what is the crush-curve of a 300 m asteroid with 50% 
porosity?), while they have a significant effect on the mo¬ 
mentum transfer efficiency. 

5.5. Selective sampling in catastrophic disruptions 

The catastrophic disruption of a large asteroid as a re¬ 
sult of a collision with a smaller projectile and the subse¬ 
quent reaccumulation of fragments as a result of their mu¬ 
tual gravitational attractions have been simulated numeri¬ 
cally during the last decades, which allowed to be repro¬ 
duced successfully the formation of asteroid families (see 
Michel et al., this volume). It is generally found that most 
large bodies formed during a catastrophic disruption consist 
of aggregates formed by reaccumulation of smaller frag¬ 
ments (e.g. Michel et al. 2001). The original location 
within the parent body of the small pieces that eventually 
reaccumulate to form the largest offspring of a disruption as 
a function of the original internal structure of the disrupted 
asteroid is interesting to determine for several reasons. If 
reaccumulation is a random process, we expect the parti¬ 
cles of a given large fragment to originate from uncorrelated 
regions within the parent body. Conversely, if the initial 
velocity field imposed by the fragmentation process deter¬ 
mines the reaccumulation phase, the particles belonging to 
the same fragment should originate from well defined areas 
in-side the parent body. In addition, the position and extent 
of these regions provide indications about the mixing occur¬ 
ring as a result of the reaccumulation process. Motivated by 
the question of the origin of ureilites, which, in some pet- 
rogenetic models, are inferred to have formed at particular 
depths within their parent body, Michel et al. (2014) started 
to investigate this problem by simulation the fragmentation 
(and reaccumulation) of a large body with a diameter fixed 
at 250 km. This study was performed using a hybrid ap¬ 
proach with a SPH shock physics code ( Benz and Asphaug 
1994,1995; Jutzi et al. 2008, Jutzi, 2014) and the n-body 
code pkdgrav as described in section 4. They considered 
four kinds of internal structures that may represent the in¬ 
ternal structure of a large body in various early stages of the 
Solar System evolution: fully molten, half molten (i.e., a 26 
km-deep outer layer of melt containing half of the mass), 
solid except a thin molten layer (8 km thick) centered at 10 
km depth, and fully solid. The properties of basalt material 
were assumed for the solid component. They focused on 
the three largest off spring that had enough reaccumulated 
pieces to consider. As found by Michel et al. (2004) who 
considered fully solid and pre-shattered bodies in the con¬ 
text of family formation, they found that the particles that 
eventually reaccumulate to form the largest reaccumulated 
bodies retain a memory of their original locations in the par¬ 
ent body. In other words, most particles in each reaccumu- 



Fig. 10. — Original positions of the material that will reaccu¬ 
mulate and form the largest, second, and third largest fragments 
(plume-shaped objects shown in different levels of gray). The re¬ 
sults for four different internal structures are shown: top left: fully 
molten; top right: half-molten; bottom left: with a thin molten 
layer (10% of the total mass) at 10 km depth (the third largest is 
hidden behind the visible two largest); bottom right: fully solid. 
The transparent gray bodies indicate the original target and im- 
pactor, respectively. The impactor moves vertically down. In the 
investigated greatly disruptive regimes, all the material that is in 
not in the largest fragments (i.e., the largest fraction of the parent 
body) is blown away, i.e., it will not reaccumulate (or only in very 
small fragments). Inspired from Michel et al. (2014). 

lated body are clustered from the same original region, even 
if their reaccumulations take place far away. However, they 
also found that the extent of the original region varies con¬ 
siderably depending on the internal structure of the parent 
(Fig. p~0] >. In particular, it seems to shrink with the solidity 
of the body. Although the covered parameter space in this 
first study was limited, this sort of investigation can give 
some constraints on the internal structure of parent bodies 
of some meteorites. 

6. CONCLUSIONS 

Numerical simulations provide an important tool that al¬ 
low us to probe regimes unreachable by experimental meth¬ 
ods. Collisions among asteroids take place in those regimes. 
To realistically model these events, the combined effects of 
gravity, strength, porosity as well as shape and structural 
properties need to be taken into account. Therefore, high 
fidelity physical models are required, and the modeling as¬ 
teroid collisions is extremely complex. 

Our knowledge of the asteroid properties which are rel¬ 
evant in terms of impact modeling is still quite limited. 
Moreover, there are known important shortcomings in the 
modeling that has been used to date, although that model¬ 
ing is improving. It is important to keep that in mind in 
assessing the meaning of any numerical simulation. 

An important step in the development of numerical 
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methods (which includes the implementation of complex 
material models) is the validation against laboratory ex¬ 
periments. Although only very briefly mentioned in this 
chapter, experimental studies greatly contributed to our un¬ 
derstanding of the collisional process at small laboratory 
scales and always give a necessary point of reference to 
test numerical methods. However, in comparisons to ex¬ 
periments it is important to include more than one single 
scalar value such as the crater size or the fragment mass 
from one experiment. In order for a numerical model to 
be a predictive tool, it should be able to reproduce multi¬ 
ple experiments in different regimes without adjusting any 
parameters. 

Since the publication of Asteroids III a decade ago, there 
have been major advances in the modeling of asteroid col¬ 
lisions and impact processes. The numerical simulations 
will undoubtedly become of higher fidelity as our models 
improve, guided by experimental, theoretical and scaling 
results. 
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